Paramétrage de l’environnement de travail

#Import des librairies utilisées
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(ggplot2)
library(tidyr)
library(tseries)
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
library(forecast)
library(TTR)

#Définition du répertoire de travail
path = "/home/utilisateur/Documents/Computer Science/ML/M1/time_series/time-series_projet/Dataset"

1. Examine your data

#Chargement du fichier day.csv et transformation en dataframe
data_day <- read.csv(file = paste(path, "/day.csv", sep=""))
df_data_day <- data.frame(data_day)

#Conversion de la variable date en format date
df_data_day$dteday <- as.Date(data_day$dteday)

#Renommage des saisons
df_data_day$season <- factor(format(df_data_day$season, format="%A"),
                          levels = c("1", "2","3","4"), 
                          labels = c("Printemps", "Eté", "Automne", "Hiver"))

#Création de nouvelles variables numériques dé-normalisées pour faciliter l'interprétation
denormalize <- function(x, minval, maxval) {
  return (x*(maxval-minval) + minval)}
df_data_day$true_temp <- Map(denormalize, df_data_day$temp, -8, 39)
df_data_day$true_temp <- df_data_day$temp*41
df_data_day$true_atemp <- df_data_day$atemp*50
df_data_day$true_windspeed <- df_data_day$windspeed*67
df_data_day$true_hum <- df_data_day$hum*100

#Visualisation et résumé statistique
head(df_data_day)
##   instant     dteday    season yr mnth holiday weekday workingday weathersit
## 1       1 2011-01-01 Printemps  0    1       0       6          0          2
## 2       2 2011-01-02 Printemps  0    1       0       0          0          2
## 3       3 2011-01-03 Printemps  0    1       0       1          1          1
## 4       4 2011-01-04 Printemps  0    1       0       2          1          1
## 5       5 2011-01-05 Printemps  0    1       0       3          1          1
## 6       6 2011-01-06 Printemps  0    1       0       4          1          1
##       temp    atemp      hum windspeed casual registered  cnt true_temp
## 1 0.344167 0.363625 0.805833 0.1604460    331        654  985 14.110847
## 2 0.363478 0.353739 0.696087 0.2485390    131        670  801 14.902598
## 3 0.196364 0.189405 0.437273 0.2483090    120       1229 1349  8.050924
## 4 0.200000 0.212122 0.590435 0.1602960    108       1454 1562  8.200000
## 5 0.226957 0.229270 0.436957 0.1869000     82       1518 1600  9.305237
## 6 0.204348 0.233209 0.518261 0.0895652     88       1518 1606  8.378268
##   true_atemp true_windspeed true_hum
## 1   18.18125      10.749882  80.5833
## 2   17.68695      16.652113  69.6087
## 3    9.47025      16.636703  43.7273
## 4   10.60610      10.739832  59.0435
## 5   11.46350      12.522300  43.6957
## 6   11.66045       6.000868  51.8261
summary(df_data_day)
##     instant          dteday                 season          yr        
##  Min.   :  1.0   Min.   :2011-01-01   Printemps:181   Min.   :0.0000  
##  1st Qu.:183.5   1st Qu.:2011-07-02   Eté      :184   1st Qu.:0.0000  
##  Median :366.0   Median :2012-01-01   Automne  :188   Median :1.0000  
##  Mean   :366.0   Mean   :2012-01-01   Hiver    :178   Mean   :0.5007  
##  3rd Qu.:548.5   3rd Qu.:2012-07-01                   3rd Qu.:1.0000  
##  Max.   :731.0   Max.   :2012-12-31                   Max.   :1.0000  
##       mnth          holiday           weekday        workingday   
##  Min.   : 1.00   Min.   :0.00000   Min.   :0.000   Min.   :0.000  
##  1st Qu.: 4.00   1st Qu.:0.00000   1st Qu.:1.000   1st Qu.:0.000  
##  Median : 7.00   Median :0.00000   Median :3.000   Median :1.000  
##  Mean   : 6.52   Mean   :0.02873   Mean   :2.997   Mean   :0.684  
##  3rd Qu.:10.00   3rd Qu.:0.00000   3rd Qu.:5.000   3rd Qu.:1.000  
##  Max.   :12.00   Max.   :1.00000   Max.   :6.000   Max.   :1.000  
##    weathersit         temp             atemp              hum        
##  Min.   :1.000   Min.   :0.05913   Min.   :0.07907   Min.   :0.0000  
##  1st Qu.:1.000   1st Qu.:0.33708   1st Qu.:0.33784   1st Qu.:0.5200  
##  Median :1.000   Median :0.49833   Median :0.48673   Median :0.6267  
##  Mean   :1.395   Mean   :0.49538   Mean   :0.47435   Mean   :0.6279  
##  3rd Qu.:2.000   3rd Qu.:0.65542   3rd Qu.:0.60860   3rd Qu.:0.7302  
##  Max.   :3.000   Max.   :0.86167   Max.   :0.84090   Max.   :0.9725  
##    windspeed           casual         registered        cnt      
##  Min.   :0.02239   Min.   :   2.0   Min.   :  20   Min.   :  22  
##  1st Qu.:0.13495   1st Qu.: 315.5   1st Qu.:2497   1st Qu.:3152  
##  Median :0.18097   Median : 713.0   Median :3662   Median :4548  
##  Mean   :0.19049   Mean   : 848.2   Mean   :3656   Mean   :4504  
##  3rd Qu.:0.23321   3rd Qu.:1096.0   3rd Qu.:4776   3rd Qu.:5956  
##  Max.   :0.50746   Max.   :3410.0   Max.   :6946   Max.   :8714  
##    true_temp        true_atemp     true_windspeed      true_hum    
##  Min.   : 2.424   Min.   : 3.953   Min.   : 1.500   Min.   : 0.00  
##  1st Qu.:13.820   1st Qu.:16.892   1st Qu.: 9.042   1st Qu.:52.00  
##  Median :20.432   Median :24.337   Median :12.125   Median :62.67  
##  Mean   :20.311   Mean   :23.718   Mean   :12.763   Mean   :62.79  
##  3rd Qu.:26.872   3rd Qu.:30.430   3rd Qu.:15.625   3rd Qu.:73.02  
##  Max.   :35.328   Max.   :42.045   Max.   :34.000   Max.   :97.25

How do the temperatures change across the seasons? What are the mean and median temperatures?

g1 <- ggplot(df_data_day, aes(x=season, y=true_temp)) + 
  geom_boxplot()
g1 + ggtitle("Boxplots des températures par saison") +
  xlab("Saison") + ylab("Température (°C)")

La période la plus chaude correspond à l’automne et la plus froide au printemps.

What are the mean and median temperatures?

mean(df_data_day$true_temp)
## [1] 20.31078
median(df_data_day$true_temp)
## [1] 20.43165

La température moyenne est de 20,3°C (0,495° normalisé) et la température médiane est de 20,4°C (0,498° normalisé).

Is there a correlation between the temp/atemp/mean.temp.atemp and the total count of bike rentals?

#Visualisation de la variable count par rapport à la variable temp
g2 <- ggplot(df_data_day, aes(true_temp, cnt)) +
  geom_point() +
  geom_smooth()
g2 + ggtitle("Locations en fonction de la température") +
  xlab("Température (°C)") + ylab("Nombre de locations")
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'

#Visualisation de la variable count par rapport à la variable atemp
g3 <- ggplot(df_data_day, aes(true_atemp, cnt)) +
  geom_point() +
  geom_smooth()
g3 + ggtitle("Locations en fonction de la température ressentie") +
  xlab("Température ressentie (°C)") + ylab("Nombre de locations")
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'

#Création de la variable mean_true_temp_atemp (moyenne quotidienne de la température et de la température ressentie)
df_data_day$mean_true_temp_atemp <- rowMeans(select(df_data_day, c(true_temp, true_atemp)))

#Visualisation de la variable count par rapport à la variable mean_true_temp_atemp
g4 <- ggplot(df_data_day, aes(mean_true_temp_atemp, cnt)) +
  geom_point() +
  geom_smooth()
g4 + ggtitle("Locations en fonction de la température moyenne \n (réelle et ressentie)") +
  xlab("Température moyenne (°C)") + ylab("Nombre de locations")
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'

D’après les nuages de points il semble exister une relation positive entre le nombre de locations et la température jusqu’à ce que les températures soient trop chaudes, au delà de 30°C environ. La relation devient alors négative.

#Calcul des coefficients de corrélation
print(c("R2 locations/température:", cor(df_data_day$true_temp, df_data_day$cnt)))
## [1] "R2 locations/température:" "0.627494009033492"
print(c("R2 locations/temp. ressentie:", cor(df_data_day$atemp, df_data_day$cnt)))
## [1] "R2 locations/temp. ressentie:" "0.631065699849181"
print(c("R2 locations/temp. moyenne:", cor(df_data_day$mean_true_temp_atemp, df_data_day$cnt)))
## [1] "R2 locations/temp. moyenne:" "0.630660733761854"

Les corrélations sont toutes d’environ 0,63. Ces variables sont donc corrélées.

What are the mean temperature, humidity, windspeed and total rentals per months?

df_data_day %>%
  group_by(data_day$mnth) %>%
  summarise_at(vars(true_temp, true_hum, true_windspeed, cnt), 
               list(mean = mean))
## # A tibble: 12 × 5
##    `data_day$mnth` true_temp_mean true_hum_mean true_windspeed_mean cnt_mean
##              <int>          <dbl>         <dbl>               <dbl>    <dbl>
##  1               1           9.69          58.6                13.8    2176.
##  2               2          12.3           56.7                14.5    2655.
##  3               3          16.0           58.8                14.9    3692.
##  4               4          19.3           58.8                15.7    4485.
##  5               5          24.4           68.9                12.3    5350.
##  6               6          28.0           57.6                12.4    5772.
##  7               7          31.0           59.8                11.1    5564.
##  8               8          29.1           63.8                11.6    5664.
##  9               9          25.3           71.5                11.1    5767.
## 10              10          19.9           69.4                11.7    5199.
## 11              11          15.1           62.5                12.3    4247.
## 12              12          13.3           66.6                11.8    3404.

Is temperature associated with bike rentals (registered vs casual)?

model_registered <- lm(registered~true_temp, data = df_data_day)
print(c("R2 for bike rentals registered:", summary(model_registered)$r.squared))
## [1] "R2 for bike rentals registered:" "0.291612923597918"
model_casual <- lm(casual~true_temp, data = df_data_day)
print(c("R2 for bike rentals casual:", summary(model_casual)$r.squared))
## [1] "R2 for bike rentals casual:" "0.295158223619129"
model_ttl <- lm(cnt~true_temp, data = df_data_day)
print(c("R2 for bike rentals total:", summary(model_ttl)$r.squared))
## [1] "R2 for bike rentals total:" "0.393748731372924"
#Affichage des résultats du modèle complet (registered et casual)
summary(model_ttl)
## 
## Call:
## lm(formula = cnt ~ true_temp, data = df_data_day)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4615.3 -1134.9  -104.4  1044.3  3737.8 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 1214.642    161.164   7.537 1.43e-13 ***
## true_temp    161.969      7.444  21.759  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1509 on 729 degrees of freedom
## Multiple R-squared:  0.3937, Adjusted R-squared:  0.3929 
## F-statistic: 473.5 on 1 and 729 DF,  p-value: < 2.2e-16
plot(model_ttl, col = "green")

Il y a un lien car la p-value est inférieure à 5% et donc l’hypothèse nulle doit être rejeté mais ce lien est faible. En effet la variable température explique moins de 40% (R2 = 0,39) des locations. D’autres variables doivent être prises en compte dans le modèle.

In the following, we you build a predictive model ff the number of bike sharing by day (daily variable names __cnt_)

Plot the cnt vs dteday and examine its patterns and irregularities

g5 <- ggplot(df_data_day, aes(dteday, cnt)) + 
  geom_point() +
  geom_smooth(method=lm)
g5 + ggtitle("Locations en fonction de la date") +
  xlab("Date") + ylab("Nombre de locations")
## `geom_smooth()` using formula = 'y ~ x'

Il y a une tendance générale avec une croissance entre début 2011 et fin 2012. Il y a également une saisonalité : on remarque une hausse des locations durant l’été et l’automne et une baisse durant l’hiver et le printemps.

Clean up any outliers or missing values if needed

#Recherche des outliers dans la variable cnt
boxplot.stats(df_data_day$cnt)$out
## integer(0)
#Recherche des valeurs NA et Null
sum(is.na(df_data_day))
## [1] 0
is.null(df_data_day)
## [1] FALSE

Il n’y a ni outliers, ni valeurs NA, ni valeurs Null

2. Now you will be using the smoothed version of cnt: choose the smoothing method and justify your choice.

Add the right frequency to your smoothed time series et justify your choices

#Création d'un time series
ts_cnt <- ts(df_data_day$cnt, start=c(2011,1), end=c(2013,1), frequency=365)

Nous travaillons avec données quotidiennes entre le 01/01/2011 (start=c(2011,1)) et le 31/12/2012 (end=c(2013,1)), le paramètre frequency doit donc être défini sur 365 (365 jours de l’année).

#Choix de la méthode de lissage
ts_cnt_HW <- HoltWinters(ts_cnt)

La moyenne mobile ne prend pas en compte les tendances et les saisonnalités. Or nos données en comporte. La méthode de lissage exponentiel de Holt-Winters a des paramètres permettant de les prendre en compte et est donc plus adaptée.

What could you tell about this new time series in term of stationarity and seasonality? Justify your conclusions.

plot(ts_cnt_HW)

plot(fitted(ts_cnt_HW))

La décomposition de la série temporelle permet d’estimer les effets de la tendance et de la saisonnalité. Sur le graphique on voit que la tendance a été retirée et la saisonalité a été atténuée.

3 Could you model the smoothed time series using ARIMA model?

#Affichage de l'ACF et de la PACF
ts_cnt_HW_fitted <- as.data.frame(fitted(ts_cnt_HW))
ts_cnt_HW_val <- ts_cnt_HW_fitted$xhat
ggtsdisplay(ts_cnt_HW_val)

La PACF décroit rapidement ce qui suggère un modèle MA(q). Mais sur l’ACF de trop nombreuses valeurs se suivent en dehors de l’intervalle de confiance. Il faut différencier la série car elle n’est pas stationaire.

#Vérification du nombre de différenciation nécessaire
ndiffs(ts_cnt_HW_val)
## [1] 1

Il faut différencier la série une fois.

#Différentiation
ts_cnt_HW_diff <- diff(ts_cnt_HW_val)
plot(ts_cnt_HW_diff)

Le bruit est réparti aléatoirement autour de 0.

#Nouvel affichage de l'ACF et de la PACF
ggtsdisplay(ts_cnt_HW_diff)

Cette fois l’ACF décroit rapidement, ce qui suggère un modèle AR(p) avec p = 5.

#Réalisation du test Augmented Dickey-Fuller afin d'obtenir des informations plus précises
adf.test(ts_cnt_HW_diff, alternative="stationary")
## Warning in adf.test(ts_cnt_HW_diff, alternative = "stationary"): p-value smaller
## than printed p-value
## 
##  Augmented Dickey-Fuller Test
## 
## data:  ts_cnt_HW_diff
## Dickey-Fuller = -10.459, Lag order = 7, p-value = 0.01
## alternative hypothesis: stationary

La p-value est inférieure à 0,05 ce qui indique que la série est stationaire. Après différenciation nous avons donc un modèle AR(5).

ts_cnt_HW_arima <- arima(ts_cnt_HW_diff, order=c(5,0,0))
ts_cnt_HW_arima
## 
## Call:
## arima(x = ts_cnt_HW_diff, order = c(5, 0, 0))
## 
## Coefficients:
##           ar1      ar2      ar3      ar4      ar5  intercept
##       -0.5264  -0.4173  -0.3613  -0.3099  -0.2114     3.0679
## s.e.   0.0511   0.0558   0.0569   0.0558   0.0512    18.7755
## 
## sigma^2 estimated as 1018331:  log likelihood = -3042.92,  aic = 6099.85

Avec 5 paramètres, ce modèle est complexe. Il faut en chercher un autre qui l’est moins.

4. Forecasting with ARIMA Models

I-Fit an ARIMA model on de-seasonal cnt (remove the season of cnt before fitting the model)

#Retrait de l'effet de la saison
ts_cnt_stl <- stl(ts_cnt,"periodic")
ts_cnt_season_adj <- seasadj(ts_cnt_stl) 
plot(ts_cnt_season_adj)

On constate qu’il n’y a plus de saisonnalité mais qu’il reste une tendance.

#Vérification du nombre de différenciation nécessaire
ndiffs(ts_cnt_season_adj)
## [1] 1

Il faut différencier la série une fois.

#Différentiation
ts_cnt_season_adj_diff <- diff(ts_cnt_season_adj)
plot(ts_cnt_season_adj_diff)

L’effet de la tendance a été retiré.

#Nouveau test Augmented Dickey-Fuller
adf.test(ts_cnt_season_adj_diff)
## Warning in adf.test(ts_cnt_season_adj_diff): p-value smaller than printed p-
## value
## 
##  Augmented Dickey-Fuller Test
## 
## data:  ts_cnt_season_adj_diff
## Dickey-Fuller = -13.998, Lag order = 8, p-value = 0.01
## alternative hypothesis: stationary

La p-value est plus petite que 0,05, la série est donc stationnaire.

II-Fit an ARIMA with Auto-ARIMA

#Auto-ARIMA
ts_cnt_auto_arima <- auto.arima(ts_cnt_season_adj)
ts_cnt_auto_arima
## Series: ts_cnt_season_adj 
## ARIMA(1,1,1) 
## 
## Coefficients:
##          ar1      ma1
##       0.2553  -0.9114
## s.e.  0.0418   0.0179
## 
## sigma^2 = 401844:  log likelihood = -5745.37
## AIC=11496.75   AICc=11496.78   BIC=11510.53

L’auto-arima a trouvé un modèle de paramètres 1,1,1

#Affichage des résidus
tsdisplay(ts_cnt_auto_arima$residuals, lag=20)

Il y a plusieurs valeurs hors de l’intervalle de confiance, donc qui ne sont pas du bruit blanc et non expliquée par le modèle.

#Test de Shapiro
shapiro.test(ts_cnt_auto_arima$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  ts_cnt_auto_arima$residuals
## W = 0.95282, p-value = 1.515e-14

En effet, le résultat du test de Shapiro indique que les résidus ne suivent pas une distribution normale (la p value est plus petite que 0,05). Il reste donc de l’information non traitée par le modèle.

#Test de Ljung-Box
Box.test(ts_cnt_auto_arima$residuals, type="Ljung-Box")
## 
##  Box-Ljung test
## 
## data:  ts_cnt_auto_arima$residuals
## X-squared = 0.0070763, df = 1, p-value = 0.933

Toutefois d’après le test de Ljung-Box, les résidus ne sont pas auto-corrélés (la p-value est plus grande que 0,05).

III-Evaluate and iterate

Le modèle auto-ARIMA ayant peu de paramètres, un modèle plus complexe pourrait expliquer d’avantage les informations restées dans le bruit. Nous allons évaluer 100 modèles en testant différentes valeurs pour les paramètres p et q afin de rechercher le modèle qui minimise le score AIC.

Compare model errors and fit criteria
#Itération sur les paramètres p et q entre 0 et 9 
aic.values <- c()
for (p in (0:9)){
  for (q in (0:9)){
    ts_cnt_season_adj_arima <- arima(ts_cnt_season_adj, order = c(p,1,q), method="ML")
    aic.values <- c(aic.values, ts_cnt_season_adj_arima$aic)    
    
  }

}

#Modèle ayant le plus petit AIC
print(c("Modèle avec le plus petit AIC:", which.min(aic.values)))
## [1] "Modèle avec le plus petit AIC:" "98"
#Valeur de son AIC
print(c("Valeur de l'AIC la plus basse:", aic.values[98]))
## [1] "Valeur de l'AIC la plus basse:" "11431.6205683934"

Le modèle qui minimise l’AIC est donc un modèle ARIMA d’ordre (9,1,7).

#Création du modèle de paramètres p = 9 et q = 8.
ts_cnt_season_adj_arima <- arima(ts_cnt_season_adj, order = c(9,1,7), method="ML")
## Warning in log(s2): NaNs produced

## Warning in log(s2): NaNs produced

## Warning in log(s2): NaNs produced
## Warning in arima(ts_cnt_season_adj, order = c(9, 1, 7), method = "ML"): possible
## convergence problem: optim gave code = 1
#Affichage du modèle
ts_cnt_season_adj_arima
## 
## Call:
## arima(x = ts_cnt_season_adj, order = c(9, 1, 7), method = "ML")
## 
## Coefficients:
##           ar1     ar2      ar3     ar4      ar5      ar6      ar7     ar8
##       -0.3747  0.1433  -1.1393  -0.446  -0.0106  -0.1185  -0.0272  0.0441
## s.e.   0.1686  0.2158   0.1184   0.181   0.2608   0.1634   0.0617  0.0621
##          ar9      ma1      ma2     ma3      ma4      ma5     ma6      ma7
##       0.1999  -0.2992  -0.5366  1.1572  -0.3980  -0.3814  0.1207  -0.2981
## s.e.  0.0492   0.1703   0.1518  0.1798   0.2069   0.1867  0.1476   0.1286
## 
## sigma^2 estimated as 348471:  log likelihood = -5698.81,  aic = 11431.62
#Affichage des résidus
tsdisplay(ts_cnt_season_adj_arima$residuals, lag=20)

On constate que si ce modèle optimisé a un AIC plus petit que le modèle auto-ARIMA (respectivement 11431.73 et 11496.75), la différence n’est pas grande et ce modèle n’explique pas mieux les résidus. En revanche ce modèle est beaucoup plus complexe que le modèle auto-ARIMA de paramètres 1,1,1. Il est donc préférable de continuer avec le modèle auto-ARIMA pour éviter l’overfitting.

Calculate forecast using the chosen model
#Prédictions sur 30 jours avec le modèle auto-ARIMA
plot(forecast(ts_cnt_auto_arima, h=30))

Plot both the original and the forecasted time series
#Original
plot(ts_cnt_season_adj, col="red")
legend(x="bottomright", legend=c("Original", "Fitted"), col=c("red", "blue"), lty=1:2, cex=0.6)

#Forecasted
lines(fitted(ts_cnt_auto_arima), col="blue")

Les prédictions sont plutôt fidèles. Les gros pics ont été lissés.

IV Forecasting

Split the data into training and test times series
end_time = time(ts_cnt_season_adj)[700]
train_set <- window(ts_cnt_season_adj, end=end_time)
test_set <- window(ts_cnt_season_adj, start=end_time)
Fit an Arima model, manually and with Auto-Arima on the training part
#Initialisation du modèle ARIMA (9,1,7)
manual_fit <- Arima(train_set, order=c(9, 1, 7))
manual_forecast <- forecast(manual_fit, h=30)

#Calcul de la précision du modèle
print(paste("Précision (RMSE) du modèle manuel: ",
            accuracy(manual_forecast, test_set)[2,"RMSE"]))
## [1] "Précision (RMSE) du modèle manuel:  602.038169309621"
#Initialisation d'un modèle avec auto-ARIMA
auto_fit <- auto.arima(train_set, seasonal=FALSE)
auto_forecast <- forecast(auto_fit, h=30)

#Calcul de la précision du modèle
print(paste("Précision (RMSE) du modèle auto-ARIMA: ", 
            accuracy(auto_forecast, test_set)[2,"RMSE"]))
## [1] "Précision (RMSE) du modèle auto-ARIMA:  694.107740668932"
#Affichage et comparaison des résultats
plot(ts_cnt_season_adj, col="black", main = "Comparaison des différents modèles")
legend(x="bottomright", legend=c("Données", "Modèle manuel", "Modèle auto"), col=c("black", "red", "green"), lty=1:2, cex=0.6)
lines(fitted(manual_fit), col="red")
lines(fitted(auto_fit), col="green")

D’après la mesure RMSE et comme l’indique le graphique, le modèle ARIMA de paramètres (9,1,7) minimise d’avantage les erreurs que le modèle manuel.

Forecast the next 25 observation and plot the original ts and the forecasted one
#Rappel de la précision du modèle manuel
print("Précision du modèle manuel:")
## [1] "Précision du modèle manuel:"
accuracy(manual_forecast)
##                    ME     RMSE     MAE       MPE     MAPE      MASE
## Training set 28.95805 598.7046 435.867 -1.568717 10.99067 0.1815598
##                      ACF1
## Training set -0.002659762
#Rappel de la précision du modèle auto-ARIMA
print("Précision du modèle auto-ARIMA:")
## [1] "Précision du modèle auto-ARIMA:"
accuracy(auto_forecast)
##                     ME     RMSE      MAE       MPE     MAPE      MASE
## Training set -2.446428 635.8652 455.2605 -2.656829 11.65282 0.1896381
##                     ACF1
## Training set 0.005113433
#Rappel de la précision du modèle Holt-Winters
print("Précision du modèle Holt-Winters:")
## [1] "Précision du modèle Holt-Winters:"
accuracy(ts_cnt_HW_arima)
##                   ME     RMSE      MAE      MPE     MAPE      MASE         ACF1
## Training set 2.77562 1009.124 740.1123 30088.77 30167.76 0.5646584 -0.009851498

La comparaison des scores nous encourage à utiliser le modèle manuel, plutôt que le modèle auto-ARIMA ou Holt-Winters. Par exemple le RMSE du modèle manuel est d’environ 599 contre 636 pour le modèle auto et 1009 pour le modèle Holt-Winters, ces derniers se trompent donc davantage.

#Prédictions sur 25 jours avec le modèle auto-ARIMA
plot(forecast(auto_forecast, h=25))

#Prédictions sur 25 jours avec le modèle manuel (9,1,7)
plot(forecast(manual_forecast, h=25))

#Test de Shapiro
shapiro.test(manual_forecast$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  manual_forecast$residuals
## W = 0.95935, p-value = 5.265e-13
#Test de Ljung-Box
Box.test(manual_forecast$residuals, type="Ljung-Box")
## 
##  Box-Ljung test
## 
## data:  manual_forecast$residuals
## X-squared = 0.0049733, df = 1, p-value = 0.9438

Le résultat du test de Shapiro indique que les résidus ne suivent pas une distribution normale (la p value est plus petite que 0,05). Toutefois d’après le test de Ljung-Box, les résidus ne sont pas auto-corrélés (la p-value est plus grande que 0,05). Comme vu précédemment il y a encore de l’information contenue dans le bruit qui n’est pas expliquée par notre modèle.

Bien que le RMSE du modèle manuel soit meilleur que celui du modèle auto-ARIMA, le modèle manuel est en revanche beaucoup plus complexe puisqu’il est de paramètres 9,1,7 contre 1,1,1. L’amélioration du RMSE apporté par tous ces paramètres supplémentaires ne vaut sans doute pas le coup par rapport à la complexité ajoutée et aux ressources de calcul nécessaires. De plus ce modèle risque plus probablement l’overfitting. Il vaudrait donc mieux utiliser le modèle auto-ARIMA.

Toutefois avec ces deux modèles, une partie de l’information reste contenue dans les résidus. Afin d’améliorer notre modèle de prédiction nous pourrions tester un modèle permettant de prendre en compte plusieurs saisonnalités (ici nous avons au moins des saisonnalités de semaine en plus des années). La fonction tbats() du package forecast le permet.

Une autre option serait d’inclure d’autres variables. Ceci peut être fait avec une régression. Une régression dynamique permettrait également de prendre en compte l’erreur et de la considérer comme une variable à expliquer et non comme un bruit comme dans une régression classique. La technique de stepwise regression pourra ensuite être utilisée pour sélectionner le modèle offrant le meilleur équilibre entre complexité et précision.